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Abstract 

This  paper  reviews  the  development  of  first-principles  based  mathematical  models  for  batteries  developed  on  a  framework  parallel  to 
computation  fluid  dynamics  (CFD),  herein  termed  computational  battery  dynamics  (CBD).  This  general-purpose  framework  makes  use  of  the 
similarity  in  the  equations  governing  different  battery  systems,  and  has  resulted  in  the  development  of  robust  models  in  a  relatively  short  time. 
Here  we  review  this  framework,  in  the  context  of  applications  to  the  coupled  modeling  of  the  thermal  and  electrochemical  behavior  of  cells, 
and  to  the  modeling  at  three  different  scales,  namely  pore-level,  cell-level  and  stack-level.  The  similarity  and  differences  of  our  approach  with 
other  research  groups  are  exemplified.  Significant  results  from  each  of  these  advanced  applications  of  modeling  are  highlighted  with  emphasis 
on  the  insights  that  can  be  gained  from  a  first-principles  model.  In  addition,  we  also  demonstrate  the  usefulness  of  a  combined  experimental- 
modeling  approach  in  describing  cells.  The  models  reviewed  here  are  expected  to  be  useful  in  predicting  the  behavior  of  advanced  batteries 
used  in  electric  vehicles  (EVs)  and  hybrid  electric  vehicles  (HEVs). 
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1.  Introduction 

The  recent  interest  in  the  development  of  environmentally 
friendly  vehicles  powered  either  fully  (electric  vehicles, 
EVs)  or  partially  (hybrid  electric  vehicles,  HEVs)  by  bat¬ 
teries  has  resulted  in  enhanced  research  into  various  aspects 
of  this  electrochemical  energy  source.  These  range  from  the 
development  of  better  designs  to  suit  the  application  (e.g. 
designs  favoring  high  power  for  HEVs  while  favoring  high 
energy  for  EVs)  to  the  development  of  new  materials 
with  superior  performance.  Parallel  to  various  experimental 
research  programs,  mathematical  models  that  describe  the 
behavior  of  batteries  and  their  interaction  with  other  devices 
in  a  vehicle  have  also  received  much  attention.  These  models 
range  from  those  that  describe  the  physics  of  the  various 
phenomena  in  the  cell  (first-principles  models  [1])  to  ones 
that  are  fit  to  experimental  data  under  various  conditions 
(equivalent  circuit  [2]  and  neural  network  models  [3]). 

While  models  that  are  trained  to  experimental  data  pro¬ 
vide  great  benefit  in  their  ability  to  fit  into  vehicle  models 
with  ease,  due  to  their  simple  construction  and  fast  compu¬ 
tational  speed,  they  possess  many  shortcomings.  Specifi¬ 
cally,  the  models  are  only  as  good  as  the  experimental  data 
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they  are  trained  to,  and  thereby  do  not  provide  the  ability 
to  extrapolate  beyond  the  range  of  this  data.  In  addition, 
changes  in  design  of  the  cell  do  not  permit  the  use  of  the 
same  models,  and  the  task  of  building  prototype  cells, 
collecting  data  and  training  the  model  has  to  be  repeated. 
More  importantly,  as  these  models  are  empirical  in  nature, 
they  provide  little,  if  any,  insight  into  the  working  of  the  cell. 
These  disadvantages  are  offset  by  the  use  of  first-principles 
models,  however  at  a  price  of  more  computational  demand 
and  added  complexity.  Therefore,  it  would  be  advanta¬ 
geous  to  develop  models  that  combine  the  strengths  of 
these  two  classes.  With  this  in  mind,  our  goal  is  to  develop 
models  based  on  the  physics  that  would  (i)  provide  insight 
into  the  operation  and  limiting  mechanisms;  (ii)  allow 
for  changes  in  design  (iii)  predict  behavior  of  cells  under 
wide  range  of  operating  conditions,  like  rates  and  tem¬ 
peratures;  (iv)  perform  these  operations  with  great  speed; 
(v)  allow  for  visualization  tools  and  user  friendly  inter¬ 
faces  and  (vi)  allow  for  integration  into  vehicle  models. 
These  six  goals  can  be  achieved  by  adapting  the  knowledge 
base  developed  over  the  last  three  decades  in  the  fluid  flow 
field  with  the  use  of  computation  fluid  dynamics  (CFD)  to 
the  field  of  batteries  (i.e.  computational  battery  dynamics, 
CBD). 

CBD  can  be  thought  to  encompass  four  key  compo¬ 
nents,  namely  (i)  physicochemical  model  development, 
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Nomenclature 

asj  interfacial  surface  area  per  unit  volume  (cm2/cm3) 
A  area  for  heat  transfer  (cm2) 
cavg  average  concentration  of  spherical  particle  = 
(1/2 F)  fgi(t)  dt(4nRi/(4/3)nRl)  (mol/cm3) 
ce  concentration  of  the  electrolyte  (mol/cm3) 
cs  surface  concentration  in  spherical  particle 
(mol/cm3) 

Cmax  maximum  concentration  (mol/cm3) 

Cp  specific  heat  (J/kg  K) 

I)  diffusion  coefficient  (cm2/s) 

F  Faraday’s  constant  (96487  C/eq.) 

h  equivalent  convective  heat  transfer  coefficient 
(W/cm2  K) 

inj  transfer  current  density  (A/cm2) 

i(t)  current  density  (A/cm2) 

IQ  initial  current  density  (A/cm2) 

k  rate  of  change  of  current  (A/cm2  s) 

/s  diffusion  length  (cm) 

q  volumetric  heat  generation  rate  (J/cm3  s) 

Rs  radius  of  the  spherical  particle  (cm) 

t  time  (s) 

T  temperature  (K) 

t/;  equilibrium  potential  (V) 

Vc  cell  volume  (cm  ) 

VCeii  cell  voltage  (V) 

Greek  letters 

KeC{  effective  solution  conductivity  (Q  1  cm-1) 

X  thermal  conductivity  (W/cm  K) 

p  density  (kg/cm3) 

creff  effective  matrix  conductivity  (Q  1  cm-1) 

(f>s  potential  in  the  matrix  phase  (V) 

</>e  potential  in  the  solution  phase  (V) 

Subscripts 

avg  average  over  cell 

amb  ambient 


(ii)  advanced  numerical  algorithms,  (iii)  material  and  prop¬ 
erty  characterization  and  (iv)  model  validation,  all  of  which 
are  illustrated  in  this  paper.  This  approach  has  proved  to  be 
very  useful  in  the  generation  of  new  models  for  various 
batteries  [4-10]  and  even  for  fuel  cell  systems  [11]  in  a  short 
time.  In  other  words,  the  general  framework  on  which  CBD 
is  based,  and  the  similar  equations  that  are  used  to  describe 
these  different  systems,  has  allowed  us  to  use  the  knowledge 
gained  in  developing  one  battery  model  in  the  other,  hence 
minimizing  the  learning  curve  considerably.  In  addition, 
graphical  user-interfaces  have  been  used  to  make  the 
models  user  friendly.  These  models  have  also  been  linked 
to  Simulink  and  have  been  successful  in  predicting  the 
behavior  of  EV  batteries  under  DST  and  SFUDS  tests  and 
compare  to  experimental  data  adequately  [12]. 


In  this  paper,  we  illustrate  CBD  by  reviewing  the  devel¬ 
opment  of  a  multi-scale  model  for  battery  systems,  which 
incorporates  a  pore-level  model  [13],  a  cell-level  model  [4-10] 
and  finally  a  stack  model  [14].  The  cell  model  developed 
here  is  a  coupled  thermal-electrochemical  model  that  not 
only  predicts  the  electrochemical  behavior  but  also  the 
thermal  excursions  in  the  cell  and  associated  changes  in 
electrochemistry.  We  first  review  the  development  of  the 
general-purpose  model  framework  that  solves  the  coupled 
thermal-electrochemical  problem  for  various  batteries.  The 
physics  included  in  the  model  and  the  subtleties/approxima¬ 
tions  incorporated  for  the  development  of  macro-homoge¬ 
neous  cell  models  are  detailed.  In  addition,  the  use  of  this 
framework  in  the  development  of  pore-level  models  using  a 
direct  numerical  simulation  (DNS)  technique  is  shown,  after 
which  the  solution  of  these  problems  using  a  general-pur¬ 
pose  solver  incorporating  advanced  solution  algorithms  is 
detailed.  Subsequently,  we  illustrate  the  cell-level  models 
using  the  Li-ion  cell  and  Ni-MH  cell  as  examples,  where  our 
combined  model/experimental  approach  is  illustrated.  We 
then  move  to  pore-level  models  and  show  the  usefulness  of 
the  DNS  technique.  Finally,  the  behavior  of  two  Ni-H2  cells 
connected  in  parallel  is  shown,  to  illustrate  our  ability  to 
simulate  stack  behavior.  Work  is  underway  to  combine  all 
three-scale  models  into  a  single  tool,  which  can  then  be  used 
to  design  advanced  batteries  at  various  levels  of  interest,  e.g. 
porous  electrodes  with  optimized  pore  structures,  single 
cells  with  certain  energy  and  power  capabilities,  and  stacks 
for  vehicular  applications  with  desirable  combinations  of 
safety,  performance,  and  life. 

2.  Mathematical  model  development 

2.1.  Macro-homogeneous  electrochemical  cell  model 

As  mentioned  in  Section  1,  the  models  developed  in  this 
paper  are  based  on  first-principles  and  incorporate  the 
physics  of  the  various  processes  occurring  in  the  cell. 
Fig.  1  illustrates  the  battery  modeled  in  this  study,  which 
consists  of  a  positive  electrode  (e.g.  nickel  hydroxide  or 
lithium  manganese  oxide  spinel),  a  negative  electrode  (e.g. 
metal-hydride  or  carbon)  and  a  separator  with  the  whole  cell 
filled  with  electrolyte  (e.g.  33%  KOH  in  Ni-MH  cells).  The 
active  material  may  be  pasted  onto  a  porous  grid,  as  in  a  Ni¬ 
MH  cell  or  on  a  planar  current  collector,  as  in  a  Li-ion  cells. 
The  electrodes  represented  in  the  figure  can  be  thought  of  as 
half  the  actual  electrode  thickness,  with  equivalent  processes 
occurring  in  the  other  half.  Symmetry  enables  us  to  neglect 
the  processes  in  the  other  half.  All  three  sections  consist  of  at 
least  two  phases,  metal  and  solution  with  a  gas  phase  in  some 
battery  systems  due  to  the  presence  of  overcharge  gassing 
reaction  (e.g.  oxygen  evolution  in  the  nickel  electrode).  The 
porous  nature  of  the  electrodes  results  in  added  ohmic  and 
diffusion  resistances  due  to  the  tortuous  ionic  and  electronic 
pathway.  In  addition,  the  existence  of  two  different  phases 
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Positive  Separator  Negative 

Electrode  Electrode 


Fig.  1.  Schematic  of  the  cell  that  is  described  in  the  cell-level  models, 
which  consists  of  two  porous  electrodes  with  a  separator  between  them. 
The  whole  cell  is  filled  with  electrolyte.  Consequently,  each  control 
volume  consists  of  matrix,  solution  and  in  some  cases  a  gas  phase.  The 
active  material  shown  in  the  inset  consists  of  many  particles  attached 
together,  one  of  which  is  simulated  in  Fig.  2. 

results  in  the  need  to  keep  track  of  the  volume  fraction  and 
the  surface  area  of  each  phase,  which  could  change  during 
charge/discharge,  due  to  changes  in  the  partial  molar  volume 
of  the  materials.  This  is  achieved  by  using  a  macro-homo¬ 
geneous  approach  where  the  detailed  geometry  of  the  pore  is 
ignored.  The  porosity  and  the  surface  area  are  averaged  over 
a  control  volume  in  the  porous  electrode  using  either  a 
porous  electrode  approach  [1,15-17]  or  a  volume  averaging 
approach  [9,18,19]  and  subsequently  used  to  correct  the 
various  balances  in  the  cell.  For  an  extensive  review  of 
flooded  porous  electrode  models,  see  De  Levie  [20],  Posey 
[21]  and  Newman  and  Tiedmann  [22].  For  more  details  on 
the  subsequent  developments  in  porous  electrode  theory,  see 
Weidner  [23]  and  Wang  et  al.  [9]. 

While  porous  electrode  theory  accounts  for  mass  transfer, 
ohmic  drops  and  reaction  in  a  volume  element,  many  battery 
electrodes  may  involve  additional  processes.  These  could 
include  solid  phase  diffusion,  liquid  phase  diffusion  in  pores, 
precipitation  and  subsequent  film  formation.  For  example,  in 
the  Li-ion  cell,  solid  phase  diffusion  of  lithium  in  the  active 
material  can  be  rate  limiting  under  certain  conditions,  hence 
requiring  models  to  describe  this  phenomenon.  This  intro¬ 
duces  an  additional  length  scale  into  the  problem,  which 
requires  a  separate  treatment.  If  the  particles  are  considered 
to  be  spherical,  then  the  diffusion  equation  needs  to  be 
solved  in  spherical  coordinates  at  every  volume  element,  in 
addition  to  solving  for  transport  in  the  v-direction  in  the 
porous  electrode.  This  can  be  accounted  for  using  a  pseudo 
two-dimensional  approach  [24-29],  or  by  numerically  inte¬ 
grating  the  time-dependent  boundary  condition  [30-32]. 
Owing  to  the  computational  tediousness  of  the  former, 
the  latter  has  been  preferred  in  the  literature.  One  such 


example  was  illustrated  by  Doyle  et  al.  [32]  where  the 
authors  numerically  solve  for  the  solid  phase  diffusion  using 
the  Duhamel’s  superposition  integral  and  incorporate  this 
into  their  porous  electrode  model. 

Although  the  exact  solution  method  [30-32]  provides 
considerable  improvement  in  speed  compared  to  the  pseudo 
two-dimensional  approach,  it  is  limited  to  restrictive 
assumptions  (e.g.  perfectly  spherical  particles  and  constant 
diffusion  coefficient  in  the  solid  phase,  where  the  exact 
integral  solution  is  possible).  In  addition,  it  results  in  an 
additional  numerical  step  at  each  control  volume,  hence 
consuming  CPU  time.  As  flexibility  and  speed  are  important 
criteria  for  our  models,  we  use  the  diffusion  length  concept, 
developed  by  Wang  et  al.  [9]  where  the  surface  and  average 
concentration  is  assumed  to  be  linearly  related  via  a  diffu¬ 
sion  length.  Mathematically  this  can  be  expressed  as  [9] 

Cs(0  =  CavgO)  +l~^  (!) 

where  the  diffusion  length,  Zs,  is  related  to  the  dimensions  of 
the  particle  and  depends  on  the  morphology  [9].  For  exam¬ 
ple,  in  spherical  particles  this  is  given  by  Rs/5.  The  term  i(t) 
represents  the  current  density  that  changes  with  time.  This 
can  be  thought  of  as  the  reaction  current  in  a  porous 
electrode,  which  changes  with  state-of-charge.  Although 
Eq.  (1)  is  physically  intuitive  and  computationally  simplis¬ 
tic,  the  assumption  that  the  surface  and  average  concentra¬ 
tions  are  linearly  dependent  on  each  other  is  valid  only  after 
the  diffusion  layer  builds  up  to  its  steady  state  value.  This 
can  be  clearly  seen  in  Fig.  2,  where  the  dimensionless 
surface  minus  the  average  concentration  is  plotted  with 
dimensionless  time  in  a  spherical  particle,  where  the  reac¬ 
tion  current  is  linearly  decreasing  with  time  with  slope  k 
(i.e.  i(t )  =  /(>  —  kt).  As  expected,  the  concentration  reaches  a 
steady  state  value  as  soon  as  the  current  is  switched  on,  after 
which  it  decreases  linearly  with  time,  mirroring  the  decrease 
in  the  current.  Also  plotted  is  the  exact  solution  obtained 
using  the  Duhamel’s  integral  method  [33],  which  shows  that 
the  concentration  starts  at  zero  and  reaches  the  steady  state 
value  after  a  time  constant,  suggesting  that  the  results 
obtained  using  Eq.  (1)  would  be  inadequate  at  short  times 
or  under  dynamic  operation  (like  pulse  or  current  interrupt). 

Considering  the  exponential  increase  in  the  concentration 
at  short  times,  which  reaches  an  asymptotic  value,  Eq.  (1) 
can  be  empirically  corrected  in  order  to  obtain  a  better 
match  to  the  Duhamel’s  solution.  An  intuitively  expressed 
correction  of  the  form 

Cs(f)  =  CavgW+^[l-e-4^]  (2) 

which  incorporates  another  time-dependent  term  in  an 
exponent  similar  to  the  time  constant  for  diffusion,  with 
a  multiplier  4/3,  has  been  seen  to  provide  good  results 
under  a  wide  range  of  operating  conditions.  Surprisingly, 
it  was  shown  by  comparing  with  the  exact  solution  that  the 
multiplier  4/3  was  valid  under  all  three  geometries  (planar. 
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Fig.  2.  Dimensionless  concentration  vs.  dimensionless  time  simulated  for  a  spherical  particle  with  symmetry  at  the  center  and  a  flux  at  the  surface  decreasing 
linearly  as  (/o  —  kt)/nF.  Generated  for  (IqRs  /  nFDcmlsf)  =  1  and  (kR:J  nFlf  cmxx)  =  0.5.  The  average  concentration  is  obtained  using  faradays  law  by  integrating 
the  decreasing  current  over  time.  The  plot  shows  the  Duhamels  solution  (solid  line)  and  comparisons  with  Eq.  (1)  (dash-dotted  line)  and  (2)  (dashed  line). 


cylindrical  and  spherical).  Fig.  2  also  shows  the  prediction 
from  Eq.  (2)  where  the  adequate  fit  to  the  Duhamel’s 
solution  is  clear.  Eq.  (2)  provides  a  simple  analytical  method 
of  adequately  describing  the  diffusion  in  the  solid  phase  with 
no  added  complexity.  However,  caution  needs  to  be  exer¬ 
cised  and  results  verified,  preferably  with  exact  analytical  or 
numerical  solutions,  before  application  to  a  system  with  new 
particle  morphologies. 

While  the  above  describes  the  method  used  to  track  the 
species  surface  concentration,  the  potential  is  estimated 
based  on  the  thermodynamics  and  kinetics  of  the  system. 
In  addition,  the  reaction  at  the  electrode  surface  results  in 
changes  in  concentration  and  generation  of  current,  which 
are  then  tracked  using  charge  balance  and  mass  balance 
equations.  As  typical  cells  contain  high  concentrations  of 
ions  (e.g.  33%  KOH  is  used  in  Ni-MH  cells)  concentrated 
solution  theory  is  used  for  this  purpose  where  in  addition  to 
accounting  for  non-idealities,  the  interaction  of  the  solute 
with  each  other  is  taken  into  account  [34-37].  Charge 
balance  is  accounted  for  using  Ohms  law,  which  is  modified 
when  used  in  the  solution  phase,  in  order  to  account  for  the 
diffusion  potential. 

2.2.  Coupled  thermal— electrochemical  model 

Incorporation  of  the  above  phenomena  results  in  the 
simulation  of  the  cell  electrochemistry  with  predictions  of 
the  concentration,  SOC,  reaction  current  and  phase  poten¬ 
tials  profiles  across  the  cell  in  addition  to  cell  voltage  and 
current.  In  order  to  generate  a  thermal  model,  these  equa¬ 


tions  are  coupled  with  an  energy  balance  at  each  control 
volume,  given  by  [8,38] 

dif>^Tl  =  V(AVr)  +  q  (3) 

at 

which  accounts  for  heat  accumulation,  conduction  and 
generation.  Assuming  a  binary  electrolyte  and  neglecting 
the  enthalpy  of  mixing  and  phase  change  effects,  the  heat 
generation  term  can  be  estimated  based  on  the  various  losses 
in  the  cell,  and  expressed  as  [6,8,38] 

_ ^  ^  Qjj. 

C1  =  "y  'as)bi;(<'/-)s  ~  <^e  —  Uj)  +  "y  'fCjinjT 

j  j 

+  (7effV<?is  V0S  +  KeffV<^e  V</>e  +  tcjfV  In  ce  V<£e  (4) 

where  the  summation  is  over  all  reactions.  The  first  term  on 
the  right  represents  the  deviation  of  the  potential  in  each 
control  volume  from  the  equilibrium  potential  (irreversible 
heat).  The  second  term  arises  from  the  entropic  effects 
(reversible  heat)  [39].  In  intercalation  electrodes,  the  rever¬ 
sible  heat  can  be  used  to  gain  insight  into  the  insertion 
mechanism  and  gauge  the  nature  of  the  insertion  sites  [41]. 
Note  that  the  first  two  terms  can  be  rewritten  as  the  deviation 
of  the  cell  potential  from  the  enthalpy  potential  or  thermo¬ 
neutral  potential  of  each  reaction  [40].  Summation  of  the 
two  terms  accounts  for  the  irreversible  and  reversible  heats 
associated  with  each  electrochemical  reaction.  The  third 
term  represents  the  ohmic  heat  arising  from  the  matrix 
phase,  while  the  last  two  terms  arise  from  ohmic  heats  in 
the  solution  phase.  Note  that  all  the  terms  that  are  needed  in 
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Eq.  (4)  are  predicted  using  an  electrochemical  model,  hence 
allowing  us  to  simulate  the  thermal-electrochemical  beha¬ 
vior  of  the  cell.  However,  as  the  temperature  of  the  cell 
changes,  the  various  controlling  parameters  in  the  cell,  like 
diffusion  coefficients  and  conductivities,  change.  These  are 
estimated  based  on  Arrhenius-type  relations  or  from  mea¬ 
sured  correlations.  Note  that  Eqs.  (3)  and  (4)  represent  the 
local  heat  generation  method  and  result  in  the  estimation  of 
the  temperature  at  every  point  in  the  cell.  When  the  cell 
temperature  is  uniform  everywhere  (small  Biot  number)  one 
can  replace  this  by  a  lumped  thermal  analysis  [8,38,42]. 
Here,  the  energy  balance  integral,  given  by 
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is  used  to  describe  the  whole  cell,  with  a  heat  generation  term 
obtained  by  integrating  Eq.  (4)  over  the  whole  cell  to  yield 


qdv 


(6) 


Only  under  the  condition  that  the  reaction  distribution  is 
uniform  across  the  porous  electrode,  and  when  no  side 
reactions  are  present,  Eq.  (6)  is  reduced  to  the  expression 
derived  by  Bernardi  et  al.  [42], 
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Therefore,  either  Eq.  (5)  is  solved  with  Eq.  (6),  wherein  the 
temperature  variation  in  the  cell  is  neglected,  or  Eq.  (5)  is 
solved  with  Eq.  (7),  whereby  the  non-uniform  reaction 
distribution  is  also  neglected.  While  the  Ni-MH  model 
neglects  the  temperature  variation,  whereby  Eqs.  (5)  and 
(6)  are  solved  [8],  the  Li-ion  model  employs  the  local  heat 
generation  method  using  Eqs.  (3)  and  (4)  [6], 


2.3.  Pore-level  model  and  direct  numerical  simulation 


Some  battery  systems,  like  the  Li-ion  cell,  have  only  5-10 
particles  across  the  thickness  in  each  electrode  [43,44], 
Under  these  conditions,  the  averaging  of  the  various  quan¬ 
tities  described  in  the  porous  electrode  approach  starts  to 
break  down.  Therefore,  a  new  kind  of  modeling  framework 
is  needed  for  such  systems.  In  addition,  as  macro-homo¬ 
geneous  models  neglect  the  detailed  morphology  of  the 
microscopic  electrode/electrolyte  interface  a  study  of  these 
microscopic  effects  will  help  in  the  development  of  better 
descriptions  of  cell  behavior  and  allow  for  interfacial  engi¬ 
neering  of  porous  electrodes.  While  a  number  of  models  exist 
that  solve  for  solution  phase  potentials  and  concentrations  in 
idealized  pore  geometries  [20,45,46],  no  first-principle 
model  for  the  current  and  concentration  distributions  in  both 
the  matrix  and  solution  phases  of  a  porous  electrode  exists. 
With  this  in  mind,  the  DNS  approach,  traditionally  used  in 
modeling  turbulence  flow  in  fluid  mechanics  and  combus¬ 
tion  in  porous  media  [47],  has  been  incorporated  into  the 
framework  of  CBD  [13]. 


Fig.  3.  Schematic  of  the  active  material  configuration  simulated  in  the 
pore-level  models  using  DNS.  Model  1  (top)  consists  of  a  continuous 
sawed  active  material,  while  model  2  (bottom)  has  a  non-reactive 
conductive  matrix  which  separates  the  active  material.  The  current 
collector  is  located  to  the  left  and  the  separator  to  the  right. 


Consider  an  idealized  geometry,  as  shown  in  Fig.  3,  con¬ 
sisting  of  matrix  and  solution  phases  in  the  carbon  electrode 
of  a  Li-ion  cell.  The  geometry  is  divided  into  a  number  of 
computational  cells  that  consist  of  either  the  matrix  or  the 
solution  phase.  This  is  different  from  porous  electrode  theory, 
where  each  control  volume  has  both  matrix  and  solution 
in  it.  Charge  and  mass  balance  expressions  are  written  for 
each  control  volume,  which  remains  the  same  irrespective  of 
whether  it  describes  the  matrix  or  solution,  with  only  the 
material  properties  varying  in  each  phase.  At  the  control 
volumes  bordering  the  metal-solution  interface,  reaction 
occurs,  which  is  accounted  for  using  Butler- Volmer  kinetics. 
The  resulting  equations,  when  solved,  provide  the  concen¬ 
tration  and  potential  distributions  in  each  phase  across  the 
electrode  with  time.  This  approach  provides  us  the  flexibility 
of  directly  simulating  complex  morphologies  of  porous 
electrodes. 

Two  different  configurations  are  illustrated,  (i)  consisting 
of  a  continuous  sawed  carbon  electrode  (Fig.  3a)  and  (ii)  the 
electrode  in  (i)  is  split  with  a  conductive  but  non-reactive 
binder  (Fig.  3b)  [13].  The  electrodes  are  immersed  in  a 
solution  of  2:1  EC:DMC.  Fig.  3  is  assumed  to  be  repre¬ 
sentative  of  one  section  of  the  porous  electrode  shown  in 
Fig.  1 ,  with  repeating  units  making  up  the  whole  electrode. 
Hence  symmetry  is  assumed  to  be  valid  at  the  top  and  bottom 
of  the  modeled  unit  system.  The  right  end  of  the  figure 
represents  the  separator,  which  is  assumed  to  be  at  a  constant 
concentration  (assumed  to  be  the  initial  concentration).  All 
model  properties  are  based  on  Doyle  et  al.  [43]. 
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3.  Model  solution 

The  electrochemical  transport  phenomena,  when  mathe¬ 
matically  described,  result  in  a  set  of  coupled  non-linear 
partial  differential  equations,  which  can  be  cast  into  a 
general  form  of  the  convective-diffusion  equation 

Q(p 

—  +  v(vd>)  =  v(rv<2>)  +  s  (8) 

at 

where  (l>  is  the  general-purpose  variable  (e.g.  concentration, 
potential),  r  the  diffusion  coefficient,  v  the  velocity  and  S  is 
the  source  term  that  incorporates  all  terms  that  do  not  fit  into 
the  other  three  terms.  The  equations  are  then  discretized 
using  the  finite  volume  method  introduced  by  Patankar  [48], 
which  results  in  a  set  of  algebraic  equations.  These  equations 
are  either  solved  iteratively  one  after  the  other  or  simulta¬ 
neously  using  a  Newton  method  coupled  with  a  GMRES 
(generalized  minimal  residual)  solver.  The  use  of  these 
advanced  methods  results  in  the  simulation  of  the  discharge 
process  of  a  two-dimensional  thermal-electrochemical 
coupled  battery  model  in  approximately  10  min. 


4.  Illustration  of  CBD:  multi-scale  modeling  of  batteries 

4.1.  Cell-level  model  with  thermal— electrochemical 
coupling 

We  begin  by  examining  the  importance  of  an  electro¬ 
chemical  model  in  predicting  the  thermal  behavior  of  cells 
by  taking  the  Li-ion  cell  as  an  example.  The  results  shown 
are  for  an  EV  cell  with  a  large  height  to  thickness  ratio 


consisting  of  a  carbon  negative  and  a  LiMn204  positive 
electrode  in  2:1  EC:DMC.  Most  of  the  cell  thermal  para¬ 
meters  were  taken  from  Baker  and  Verbrugge  [49],  while  the 
thickness  of  the  electrodes  and  the  electrochemical  para¬ 
meters  (diffusion  coefficients,  equilibrium  potentials)  were 
taken  from  Doyle  et  al.  [43].  The  temperature  dependence 
of  the  various  parameters  (e.g.  exchange  current  densities) 
were  taken  from  Botte  et  al.  [50],  More  details  on  the 
parameters  used  can  be  found  in  [6].  The  entropic  term 
(see  Eq.  (4))  has  previously  been  neglected  due  to  lack  of 
data.  However,  recently  reported  data  for  the  two  electrodes 
used  in  this  study  has  permitted  us  to  use  this  information  in 
the  present  model.  For  the  manganese  oxide  spinel,  the  data 
present  by  Thomas  et  al.  [41]  was  used,  while  for  the  carbon 
electrode  the  data  presented  by  Al  Hallaj  et  al.  [51]  as  a 
function  of  equilibrium  potential  was  used  to  estimate  it  as  a 
function  of  amount  of  lithium  intercalated. 

Consider  a  case  when  the  cell  simulated  is  situated  in  the 
middle  of  a  cell  stack.  Under  these  circumstances  no  heat 
dissipation  occurs  from  the  sides  of  the  cell  and  all  dissipa¬ 
tion  occurs  from  the  top,  through  the  tabs.  Fig.  4a,  which 
shows  the  temperature  distribution  during  a  2  C  discharge  a 
2.9  Ah  Li-ion  cell  at  50%  SOC,  suggests  that  the  large  aspect 
ratio  can  lead  to  significant  two-dimensional  effects.  As  the 
top  of  the  cell  is  exposed  to  the  ambient,  the  temperature  at 
the  cell  top  can  be  considerably  lower  than  at  the  cell  bottom. 

This  distribution  of  temperature  in  the  cell  results  in  the 
reaction  current  being  different  at  different  points  in  the 
cell,  as  seen  in  Fig.  4b  [6],  Note  that  the  separator  has  zero 
reaction  current  as  no  charge  transfer  occurs  in  this  region, 
while  the  current  have  different  signs  in  the  two  electrodes. 
While  even  under  isothermal  conditions  changes  in  the 
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Fig.  4.  Temperature  (left)  in  K  and  reaction  current  (right)  in  A/cm3  contours  generated  in  a  Li-ion  cell  during  2  C  discharge  at  50%  SOC.  The  features 
incorporated  in  the  model  are  detailed  in  the  text.  See  [6]  for  more  details.  Heat  dissipation  is  assumed  to  occur  only  from  the  top  of  the  cell. 
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reaction  current  across  the  porous  electrode  could  occur 
depending  on  the  ratio  of  the  matrix  to  solution  phase 
conductivities,  the  kinetic  resistance  and  the  slope  of  the 
equilibrium  potential  with  SOC  [52],  Fig.  4b  has  the  added 
effect  of  having  different  behavior  due  to  differences  in 
temperature.  This  is  seen  more  clearly  in  the  carbon  elec¬ 
trode  where  the  changing  reaction  current  with  cell  height, 
especially  at  the  top  of  the  cell,  is  apparent.  Note  that  the 
temperature  variation  seen  in  Fig.  4a  can  decrease  consider¬ 
ably  by  decreasing  the  cell  height  or  by  choosing  current 
collectors  (which  have  two  orders  of  magnitude  larger 
thermal  conductivity)  of  larger  thickness. 

The  change  in  the  electrochemical  behavior  with  change 
in  temperature,  seen  in  Fig.  4b,  is  clearer  in  Fig.  5  where  the 
cell  voltage  is  plotted  at  different  rates  under  isothermal  and 
non-isothermal  conditions.  The  non-isothermal  case  was 
simulated  using  infinite  heat  dissipation  from  the  top  of 
the  cell  but  no  heat  dissipation  from  the  sides.  At  low  rates 
(0.01  C)  the  cell  is  practically  at  equilibrium  and  hence  heat 
generation  occurs  only  due  to  the  reversible  heat  effects. 
Flowever,  as  the  time  of  discharge  under  these  conditions  is 
large,  heat  dissipation  from  the  cell  results  in  negligible 
increase  in  the  temperature  (not  shown)  even  under  the  non- 
isothermal  case.  As  the  temperature  rise  is  small,  no  effect  is 
seen  in  the  electrochemical  behavior.  However,  as  the  rate 
increases,  the  temperature  in  the  cell  starts  to  increase.  This 
increase  results  in  a  decrease  in  the  kinetic,  mass  transfer  and 
ohmic  resistances  in  the  cell,  which  results  in  the  higher 
voltage  at  the  same  SOC.  In  addition,  the  increase  in  the 
solid  phase  diffusion  results  in  an  increase  in  the  utilization 
of  active  materials.  For  example,  during  a  2  C  discharge,  as 


much  as  25%  increase  in  utilization  is  seen  due  to  the  higher 
temperature.  In  summary,  Figs.  4  and  5  assert  that  while  the 
thermal  behavior  of  the  cell  is  dictated  by  the  electrochem¬ 
istry,  the  electrochemical  behavior  is  considerably  affected 
by  the  thermal  excursions.  The  coupling  of  the  two  gives  rise 
to  thermal  runaway. 

Considering  the  complexity  involved  in  developing  an 
electrochemical  model  for  cell,  it  is  tempting  to  avoid  this 
problem  by  using  a  lumped  thermal  model  (see  Eq.  (5))  and 
estimating  the  heat  generation  rate  from  experimental  data 
instead  of  Eq.  (6)  or  7.  For  example,  the  heat  generation  can 
be  estimated  by  either  performing  experiments  at  isothermal 
conditions  and  using  the  difference  in  the  cell  potential 
from  the  thermo-neutral  potential  [53],  or  by  performing 
experiments  using  an  isothermal  calorimeter.  In  other  words, 
this  approach  assumes  that  the  heat  generation  under  iso¬ 
thermal  conditions  is  same  as  that  under  non-isothermal 
conditions.  However,  as  seen  in  Fig.  5,  as  the  temperature  of 
the  cell  increases,  the  losses  in  the  cell  decrease,  which  in 
turn  will  decrease  the  heat  generation  rate,  hence  making 
this  assumption  suspect  [54]. 

This  can  be  clearly  seen  in  Fig.  6,  where  the  heat  energy 
from  the  cell  during  a  complete  discharge  is  plotted  with 
C-rate  for  isothermal  and  non-isothermal  cases.  The  figure 
was  generated  by  performing  discharge  simulations  at  var¬ 
ious  rates  and  integrating  the  heat  versus  time  curve  to 
calculate  the  total  thermal  energy  generated.  Fig.  6  shows 
two  distinctive  features,  namely  (i)  a  peak  in  the  energy  as  the 
discharge  rate  is  increased  and  (ii)  a  significant  difference  in 
the  energy  between  the  isothermal  and  non-isothermal  cases. 
At  low  rates,  the  cell  is  operating  at  equilibrium  conditions 


Fig.  5.  Voltage  vs.  DOD  at  different  rates  for  a  Li-ion  cell.  Graphs  are  generated  under  isothermal  conditions  (dashed  line)  and  when  heat  dissipation  occurs 
from  only  the  top  of  the  cell  (solid  line).  See  [6]  for  details  of  the  parameters  used. 
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Fig.  6.  Energy  generated  from  a  Li-ion  cell  during  discharge  at  different  C-rates.  The  graph  was  obtained  by  integrating  the  heat  generation  curve  over  the 
whole  discharge.  The  graph  compares  the  isothermal  to  the  non-isothermal  case. 


whereby  the  energy  dissipated  is  due  to  the  entropic  heat  and 
hence  a  constant.  As  the  rate  increases,  kinetic,  ohmic  and 
mass  transfer  losses  start  to  increase,  leading  to  an  increase 
in  the  energy.  However,  at  very  high  rates,  although  the  heat 
generation  is  large,  the  utilization  starts  to  decrease,  hence 
decreasing  the  time  for  heat  generation.  This  results  in  a 
decrease  in  the  total  energy  generated,  hence,  leading  to  a 
peak. 

Fig.  6  clearly  shows  that  the  non-isothermal  and  isother¬ 
mal  heat  generation  data  are  not  identical  and  can  vary  by  as 
much  as  50%  from  each  other,  whereby  predictions  of  the 
temperature  of  the  cell  would  be  inaccurate.  In  addition, 
while  the  temperature  would  be  under-predicted  at  certain 
rates,  it  would  be  over-predicted  in  others.  It  should  be  noted 
that  curves  similar  to  those  shown  in  Fig.  6  could  be 
generated  under  different  conditions  (e.g.  natural  convec¬ 
tion,  forced  convection,  adiabatic)  with  no  two  being  the 
same.  Therefore,  experimental  data  under  any  one  condition 
cannot  be  used  to  predict  thermal  behavior  under  other 
conditions,  asserting  the  need  for  a  combined  thermal- 
electrochemical  modeling  approach. 

Having  established  the  need  for  a  combined  thermal- 
electrochemical  approach,  we  now  turn  to  examples  of 
property  characterization  and  model  validation.  This  is 
illustrated  for  a  Ni-MH  cell,  where  model  to  isothermal 
experimental  comparison  was  conducted  on  a  commercial 
95  Ah,  12  V  Ni-MH  EV  battery.  The  model,  which  uses  a 
lumped  thermal  analysis,  includes  the  side  reactions  due  to 
oxygen  evolution  in  the  nickel  electrode  and  oxygen  reduc¬ 
tion  and  hydrogen  evolution  in  the  MH  electrode  [5,8]. 
Scanning  electron  microscope  images  and  BET  surface  area 


measurements  on  the  two  electrodes  showed  that  the  particle 
sizes  of  the  active  materials  (500  nm  in  the  MH  and  50  nm  in 
the  Ni  electrode)  were  small  enough  that  solid-state  diffu¬ 
sion  limitations  could  be  neglected.  In  addition,  the  cells 
were  tom  down  to  measure  the  thickness  and  porosities  of  the 
electrodes  and  separator,  which  were  then  used  in  the  model. 
Experiments  at  low  rates  (C/20)  were  assumed  to  represent 
equilibrium  conditions  and  was  used  to  estimate  the  equili¬ 
brium  potential  as  a  function  of  SOC.  It  was  seen  that  the 
nature  of  the  curve  was  very  different  from  a  Nemstian  or 
other  traditional  activity  coefficient  corrected  thermodynamic 
models  [55-57],  and  hence  was  fit  to  an  empirical  equation. 

Experiments  were  then  performed  at  a  rate  of  4  C,  where 
the  cell  showed  considerable  limitation,  although  with  little 
loss  in  capacity  at  lowered  cut-off  potentials,  asserting  to  the 
validity  of  the  assumption  that  solid  phase  diffusion  was  not 
a  limitation  mechanism  in  these  cells.  With  the  assumption 
that  the  ohmic  drops  in  the  porous  electrodes  and  separator 
were  known,  the  only  unknown  parameter  in  the  cell  is  the 
kinetics.  Therefore,  the  i0  and  a  of  the  Ni  electrode  were 
changed  in  order  to  obtain  a  adequate  match  to  the  experi¬ 
mental  data.  Using  these  values  of  the  kinetic  parameters,  the 
model  was  used  to  predict  the  behavior  of  the  cell  at  other 
rates,  as  shown  in  Fig.  7,  where  the  excellent  fit  of  the  model 
to  the  data  is  seen  at  all  four  rates,  suggesting  that  the 
phenomena  included  in  the  model  were  adequate  in  predicting 
cell  behavior  under  a  wide  range  of  operating  conditions. 

The  advantage  of  having  a  mathematical  model  that 
predicts  the  behavior  of  the  cell  is  the  ability  to  extrapolate 
to  different  operating  conditions.  One  example  of  this  is 
shown  in  Fig.  8,  where  the  cell  voltage  versus  time  is  shown 
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Fig.  7.  Experimental  vs.  model  predictions  of  discharge  curves  at  various  rates  for  a  Ni-MH  cell.  The  model  was  fit  to  the  data  generated  at  0.05  and  4  C  and 
then  used  to  predict  the  voltage  at  the  rates  shown. 


for  a  1.2  V  Ni-MH  cell  during  constant  current  charge  under 
different  heat  transfer  coefficients  ranging  from  isothermal 
to  adiabatic  conditions.  The  plot  also  shows  the  temperature 
rise  during  the  charge.  As  the  cell  shifts  from  operating  from 
isothermal  to  adiabatic  conditions,  as  expected,  the  tem¬ 
perature  of  the  cell  increases.  Under  adiabatic  conditions, 
temperatures  as  high  as  80  °C  are  predicted,  asserting  to 
the  need  for  thermal  management  systems  in  these  cells. 
Notable  in  Fig.  8  is  the  voltage  rollover  seen  close  to  the  end 
of  charge,  which  becomes  steeper  as  the  temperature  of  the 
cell  increases.  This  rollover  feature,  which  has  been  seen 
experimentally,  is  caused  by  the  oxygen  reduction  reaction 
at  the  MH  electrode.  As  a  Ni-MH  cell  is  charged,  oxygen 
evolves  from  the  positive  electrode,  according  to 

40H~  ->  02  +  2H20  +  4e~  (9) 

which  is  then  transported  to  the  negative  via  the  separator.  At 
the  negative,  the  large  driving  force  for  the  reaction  results  in 
its  consumption  by  the  reverse  of  reaction  (9).  Therefore, 
during  overcharge,  two  reactions  occur  on  the  negative, 
namely,  MH  charging  and  oxygen  reduction.  The  potential 
of  the  system  is  therefore  a  mixed  potential,  which  depends 
on  the  fraction  of  the  current  going  to  each  reaction.  When 
the  cell  temperature  increases  (i.e.  isothermal  to  adiabatic 
charge),  the  i0  of  the  oxygen  reduction  reaction  increases, 
which  in  turn  results  in  the  overpotential  for  this  reaction 
decreasing.  This  decrease  in  the  overpotential  results  in  the 
potential  of  the  negative  electrode  shifting  to  more  positive 
values  thereby  decreasing  the  cell  potential.  This  increase  in 
reaction  kinetics  is  the  cause  for  the  increased  rollover 
feature,  seen  in  Fig.  8,  as  the  temperature  increases. 


4.2.  Pore-level  model 

While  the  macro-homogeneous  models  provide  insights 
to  the  behavior  of  the  cell  from  a  macroscopic  scale,  details 
on  the  phenomena  on  a  microscopic  scale  are  lost.  In  this 
section,  we  describe  the  pore-level  models  developed  in 
order  to  capture  this  effect.  We  illustrate  this  using  the  two 
different  configurations  of  active  material,  shown  in  Fig.  3. 
We  simulate  discharge  behavior  of  a  carbon  electrode  where 
lithium  is  de-intercalated  from  the  active  material.  As  model 
geometry  2  has  non-reacting  binder  in  addition  to  active 
material,  the  total  amount  of  dischargable  material  is 
reduced.  This  is  seen  from  the  concentration  profiles  of 
Li  in  the  solid  phase  20  min  into  discharge,  shown  in  Fig.  9. 
The  plot  was  generated  at  1  C  rate  corresponding  to  model  1, 
which  is  approximately  2  C  for  model  2.  As  the  volume  of 
active  material  is  less  in  model  2  compared  to  model  1,  the 
moles  of  active  material  is  also  less.  Hence,  after  discharge 
at  the  same  current  for  the  same  time,  the  concentration  in 
model  2  is  lower  than  that  in  model  1,  as  seen  in  the  figure. 
Although  considerable  concentration  variations  are  seen  in 
the  solid  phase  concentrations  in  the  figure,  no  such  differ¬ 
ences  are  seen  in  the  liquid  phase  concentration,  indicating 
that  under  these  conditions,  liquid  phase  limitations  are 
negligible.  In  addition,  the  differences  in  the  configuration 
of  the  active  material  result  in  differences  in  the  diffusion 
flux  into  the  active  material  between  the  two  models.  Note 
that  in  these  simulations  the  binder  and  active  material  were 
assumed  to  be  the  same  electronic  conductivity. 

The  approach  shown  here  is  an  ideal  procedure  for 
evaluating  a  number  of  electrode  configurations  that  are 
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Fig.  8.  Temperature  (top)  and  voltage  (bottom)  vs.  charge  input  for  a  Ni-MH  cell  under  conditions  ranging  from  isothermal  to  adiabatic.  The  graph  shows  the 
voltage  rollover  occurring  from  oxygen  reduction  in  the  MH  electrode  [8]. 


impossible  to  describe  using  a  macro-homogeneous  model. 
These  include  the  ability  to  study  effect  of  adding  conduct¬ 
ing  binders  in  less  conducting  active  material,  irregularly 
shaped  particles  in  electrodes,  studying  particle  size  dis¬ 
tributions  as  opposed  to  a  bimodal  distribution  [58,59],  and 
studying  the  effect  of  having  graded  porosity/pore  size  in  the 
electrode  [60].  The  DNS  technique  provides  an  excellent 
method  of  evaluating  and  testing  various  electrode  config¬ 
urations  prior  to  electrode  manufacture. 

4.3.  Stack-level  model 

We  finally  examine  the  third  level  of  modeling,  namely 
stack  modeling  of  batteries.  Specifically,  we  are  interested  in 
simulating  the  behavior  of  cells  that  are  connected  with  each 


other  in  a  stack  and  the  SOC  balancing  issues  that  arise  due 
to  thermal  excursions.  Consider  two  Ni-H2  cells  connected 
together  in  parallel,  with  each  cell  at  a  different  temperature 
and  charged  at  a  constant  current  to  the  stack  [14].  This  can 
be  thought  of  as  two  cells  at  different  positions  in  a  cell 
stack.  As  the  temperatures  in  the  cells  are  different,  the 
various  resistances  in  the  cell  are  also  different,  thereby 
resulting  in  different  characteristics.  When  the  cells  are 
connected  externally  in  parallel,  the  voltages  in  the  two 
are  forced  to  be  a  constant,  which  then  results  in  differences 
in  the  current  to  each  cell,  with  their  sum  being  a  constant. 
This  is  seen  in  Fig.  10,  where  the  total  reaction  current  in  the 
two  cells  is  plotted  as  the  batteries  are  charged  [14].  The 
graph  was  made  by  using  two  cell-level  thermal-electro¬ 
chemical  models  and  imposing  the  conditions  that  the 
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Fig.  9.  Concentration  (mol/cm3)  contours  for  lithium  in  the  matrix  and  solution  during  discharge  of  a  carbon  electrode  in  a  Li-ion  cell  schematized  in  Fig.  3. 
While  the  concentration  in  solution  refers  to  Li+,  in  the  matrix  it  refers  to  Li.  The  top  frame  refers  to  model  1  and  the  bottom  frame  to  model  2. 


voltage  of  the  two  cells  remain  the  same  and  that  the  currents 
sum  up  to  the  external  current  to  the  stack. 

During  charge,  Ni  oxidation  is  the  primary  reaction  in  the 
positive  electrode,  which  on  overcharges  transitions  to  the 
oxygen  evolution  reaction.  At  low  states  of  charge,  the  cell 
with  the  higher  temperature  has  a  lower  voltage  closer  to  its 
equilibrium  value,  whereby  a  greater  current  is  passed  to  the 
cell  in  order  to  make  its  voltage  match  that  of  the  cooler  cell. 
As  charge  continues,  the  hotter  cell  reaches  full  charge  faster 
and  hence  the  current  to  the  cell  decreases  compared  to  the 
colder  cell.  However,  on  overcharge,  the  oxygen  evolution 
reaction  starts  to  occur,  which  is  more  energetic  for  the  hotter 
cell  compared  to  the  cooler  cell  hence  resulting  in  a  similar 
profile  as  that  seen  in  the  first  part  of  the  charge.  At  longer 
times,  the  Ni  reaction  reaches  completion  in  both  the  cells  and 


only  oxygen  evolution  occurs,  which  results  in  the  current 
reaching  a  steady  state.  More  details  can  be  found  in  [14]. 

Fig.  10  is  an  example  of  the  complex  interactions  that  can 
occur  in  a  cell  stack  due  to  the  cross  talk  between  the 
different  cells.  In  other  words,  simulation  of  stack  behavior, 
where  interaction  between  the  cells  is  not  considered,  would 
be  inaccurate.  In  addition.  Fig.  10  was  simulated  by  holding 
the  temperature  of  each  battery  the  same  through  the  charge; 
a  condition  that  is  not  expected  in  reality.  Efforts  are  now 
underway  to  study  cell  stacks  where  in  addition  to  electro¬ 
chemical  variations,  temperature  variations  are  also  being 
included.  A  parallel  computer  based  on  PC  clusters  is  used 
for  these  calculations.  This  study,  being  conducted  for  Li-ion 
cells,  is  expected  to  generate  useful  information  for  devel¬ 
opment  of  thermal  management  systems  for  cell  stacks. 


Fig.  10.  Reaction  current  vs.  charge  input  for  two  Ni-F^  cells,  connected  in  parallel,  one  at  283  K  and  other  at  303  K.  The  graph  is  simulated  during  charge 
when  a  constant  current  is  passed  to  the  stack  [14]. 
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5.  Conclusions 

This  paper  reviews  the  development  of  a  first-principle 
based  mathematical  models  to  describe  batteries  on  a  frame¬ 
work  parallel  to  CFD,  termed  CBD.  Specifically  two  exam¬ 
ples  of  this  approach  are  illustrated,  namely,  development  of 
thermal-electrochemical  coupled  models  and  the  develop¬ 
ment  of  multi-scale  models.  The  model  equations  and  the 
significant  results  from  a  pore-level,  cell-level  and  stack- 
level  model  are  reviewed. 

The  cell-level  models  are  used  to  show  the  importance  of 
a  combined  thermal-electrochemical  modeling  approach  to 
describe  thermal  behavior  of  batteries. 

Specifically,  experimental  data  under  one  condition  (e.g. 
isothermal)  cannot  be  used  to  predict  thermal  excursions 
under  other  conditions;  a  problem  avoided  by  the  use  of  an 
electrochemical  model.  In  addition,  the  combined  experi¬ 
mental/model  approach  is  used  to  show  the  ability  of  the 
models  to  describe  behavior  under  different  rates  using  a  Ni¬ 
MH  cell  as  example.  The  pore-level  models  were  used  to 
illustrate  the  ability  to  simulate  complex  interfacial  geome¬ 
tries,  which  are  ignored  in  cell-level  models.  Finally,  the 
ability  to  simulate  cell  stacks  was  illustrated  by  taking  two 
Ni-H2  cells  at  different  temperatures,  connected  in  parallel 
and  charged  using  as  a  constant  current,  as  an  example.  The 
importance  of  incorporating  the  interaction  between  the  two 
batteries  in  stack  simulations  is  shown. 
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